#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
######################################################################################
优化势函数: 从eamin文件读取势参数, 更新参数, 调用eamgen生成势文件, 调用lammps计算
python potfit.py --method=int --dstep=float --dmax=float --acc=float --maxfun=int --xscale=0/1 --yscale=0/1 --nproc=int --nebimage=int
Author:  Xiao-Chun Li (xcli@ipp.ac.cn)
Last Update: 2025-06-16
######################################################################################
"""

import os, sys, getopt, time, copy, platform
import numpy as np
from scipy import optimize

# 确保导入eamgen模块
try:
    if os.path.isfile("eamgen.py"):
        sys.path.append(os.getcwd())
        from eamgen import eamgen
    else:
        print("警告: eamgen.py 文件不存在")
except Exception as e:
    print(f"导入eamgen模块时出错: {e}")


######################################################################################
# 判断现有的LAMMPS路径
######################################################################################
def get_lammps_path():
    npimg = int(nproc / nebimage)
    npneb = nebimage * npimg
    lammps = ["", "", "", 1, 1]  # [0]in*运行命令, [1]RDM*运行命令, [2]NEB*运行命令, [3]RDM并行计算开关=1/0, [4]NEB并行计算开关=1/0
    if os.path.isfile("C:/Program Files/Microsoft MPI/Bin/mpiexec.exe") and os.path.isfile(
            "C:/SOFT/LAMMPS-64bit/bin/lmp.exe"):  # lammps路径@Windows
        lammps[0] = f'"C:/SOFT/LAMMPS-64bit/bin/lmp.exe" -pk omp {nproc} -sf omp'
        lammps[
            1] = f'"C:/Program Files/Microsoft MPI/Bin/mpiexec.exe" -np {nproc} C:/SOFT/LAMMPS-64bit/bin/lmp.exe -partition {nproc}x1'
        lammps[
            2] = f'"C:/Program Files/Microsoft MPI/Bin/mpiexec.exe" -np {npneb} C:/SOFT/LAMMPS-64bit/bin/lmp.exe -partition {nebimage}x{npimg}'
    elif os.path.isfile("C:/SOFT/LAMMPS-64bit/bin/lmp.exe"):  # lammps路径@Windows: 无mpi, set lammps[3] = 0, lammps[4] = 0
        lammps[0] = f'"C:/SOFT/LAMMPS-64bit/bin/lmp.exe" -pk omp {nproc} -sf omp'
        lammps[1] = "=" * 50 + f"> Warning: No mpiexe, RDM jobs run in serial module"
        lammps[2] = "=" * 50 + f"> Warning: No mpiexe, NEB jobs will pass"
        lammps[3] = 0
        lammps[4] = 0
    elif os.path.isfile("/usr/bin/lmp"):  # lammps路径@Ubuntu
        lammps[0] = f"mpirun -np {nproc} /usr/bin/lmp"
        lammps[1] = f"mpirun -np {nproc} /usr/bin/lmp -partition {nproc}x1"
        lammps[2] = f"mpirun -np {npneb} /usr/bin/lmp -partition {nebimage}x{npimg}"
    elif os.path.isfile("/public/software/apps/lammps/intelmpi/29Oct20/lmp_cpu"):  # lammps路径@合肥超算
        lammps[0] = f"mpirun -bootstrap slurm -n {nproc} /public/software/apps/lammps/intelmpi/29Oct20/lmp_cpu"
        lammps[
            1] = f"mpirun -bootstrap slurm -n {nproc} /public/software/apps/lammps/intelmpi/29Oct20/lmp_cpu -partition {nproc}x1"
        lammps[
            2] = f"mpirun -bootstrap slurm -n {npneb} /public/software/apps/lammps/intelmpi/29Oct20/lmp_cpu -partition {nebimage}x{npimg}"
    elif os.path.isfile("/public/software/lammps/lammps-2018-08-22-neb/lmp_icc_openmpi"):  # lammps路径@服务器
        lammps[
            0] = f"/public/software/mpi/openmpi/1.6.5/intel/bin/mpirun -np {nproc}  /public/software/lammps/lammps-2018-08-22-neb/lmp_icc_openmpi"
        lammps[
            1] = f"/public/software/mpi/openmpi/1.6.5/intel/bin/mpirun -np {nproc}  /public/software/lammps/lammps-2018-08-22-neb/lmp_icc_openmpi -partition {nproc}x1"
        lammps[
            2] = f"/public/software/mpi/openmpi/1.6.5/intel/bin/mpirun -np {npneb}  /public/software/lammps/lammps-2018-08-22-neb/lmp_icc_openmpi -partition {nebimage}x{npimg}"
    else:
        print("Unknow platform, please modify lammps path @ def get_lammps_path():")
        sys.exit(0)
    if npimg < 1:
        print("=" * 50 + f"> Warning: nproc={nproc}, nebimage={nebimage}, npimg={npimg}, NEB jobs will pass")
        lammps[4] = 0
    return lammps


######################################################################################
# 拟合参数获取: 如果存在control文件, 先从control文件读取默认值, 命令行参数可以覆盖control文件设置
######################################################################################
def get_fit_options():
    ##########################################################################
    # 默认参数设置
    ##########################################################################
    global xscale  # 先声明为全局变量
    
    method = 1  # 拟合方法选择, 0表示不拟合只计算, 1-4选择不同方法
    dstep = 1.0e-6  # 修改步长
    dmax = 1.0e+2  # xfit可以修改的最大值
    acc = 1.0e-6  # 收敛精度, ftol  xtol  gtol  epsfcn
    maxfun = 10000  # 最多拟合次数
    xscale = 1  # 输入参数是否归一化处理, 初始参数大于0.1的, 按1处理, 否则不变
    yscale = 0  # 计算性质是否归一化处理, 参考值<0.05时, 给出差值*500
    nproc = 8  # LAMMPS计算采用核数
    nebimage = 7  # NEB计算的镜像数, 需与NEB计算文件保持一致

    ##########################################################################
    # 如果存在control文件, 先从control文件读取默认值
    ##########################################################################
    local_vars = {}
    if os.path.isfile("control"):  # exec(object[, globals[, locals]])
        exec(open("control", "r", encoding="utf-8").read(), globals(), local_vars)
    elif os.path.isfile("control.py"):  # exec(object[, globals[, locals]])
        exec(open("control.py", "r", encoding="utf-8").read(), globals(), local_vars)
    if "method" in local_vars.keys():  method = int(local_vars["method"])
    if "dstep" in local_vars.keys():  dstep = local_vars["dstep"]
    if "dmax" in local_vars.keys():  dmax = local_vars["dmax"]
    if "acc" in local_vars.keys():  acc = local_vars["acc"]
    if "maxfun" in local_vars.keys():  maxfun = int(local_vars["maxfun"])
    if "xscale" in local_vars.keys():  xscale = int(local_vars["xscale"])
    if "yscale" in local_vars.keys():  yscale = int(local_vars["yscale"])
    if "nproc" in local_vars.keys():  nproc = int(local_vars["nproc"])
    if "nebimage" in local_vars.keys():  nebimage = int(local_vars["nebimage"])

    ##########################################################################
    # 命令行获取参数：命令行参数可以覆盖control文件设置
    ##########################################################################
    try:
        argvs, unknowargs = getopt.getopt(sys.argv[1:], shortopts="",
                                          longopts=["method=", "dstep=", "dmax=", "acc=", "maxfun=", "xscale=",
                                                    "yscale=", "nproc="])
    except getopt.GetoptError:
        print("=" * 50 + f"> ERROR: getopt {sys.argv[1:]}")
        sys.exit()
    for name, value in argvs:
        if name in ("--method"):   method = int(eval(value))
        if name in ("--dstep"):    dstep = eval(value)
        if name in ("--dmax"):     dmax = eval(value)
        if name in ("--acc"):      acc = eval(value)
        if name in ("--maxfun"):   maxfun = int(eval(value))
        if name in ("--xscale"):   xscale = int(eval(value))
        if name in ("--yscale"):   yscale = int(eval(value))
        if name in ("--nproc"):    nproc = int(eval(value))
    if len(unknowargs) > 0: print(f"Unkown argvs: {unknowargs}")
    return method, dstep, dmax, acc, maxfun, xscale, yscale, nproc, nebimage


######################################################################################
# 读取ref文件
######################################################################################
def read_ref(filename="ref"):
    if not os.path.isfile(filename):
        print("=" * 50 + f"> ERROR: file [{filename}] NOT exist.")
        sys.exit()
    ref_dict = {}  # ref_dict["Key"][0]: 参考性质, ref_dict["Key"][1]: 权重
    Preset_dict = {}  # 预设性质 Preset_dict["Key"]
    ref_Eb = []  # Eb计算公式: Eb0 = Ef1 + Ef2 - Ef3
    lines = open(filename, "r", encoding="utf-8").readlines()
    for line in lines:
        words = line.split()
        prop = words[0]
        if len(words) >= 3:
            if prop not in ref_dict.keys():
                ref_dict[prop] = [float(words[1]), float(words[2])]  # [参考性质, 权重]
            else:
                print("=" * 50 + f"> ERROR: {prop} Repeat Reference data")
                sys.exit()
        if "PresetValue" in line:
            if prop not in Preset_dict.keys():
                Preset_dict[prop] = float(words[1])
            else:
                print("=" * 50 + f"> ERROR: Repeat PresetValue: {prop}")
                sys.exit()
        if len(words) >= 6 and words[3] in ["+-", "-"]: ref_Eb.append(words)
    return ref_dict, ref_Eb, Preset_dict


######################################################################################
# 从eamin读取初始势参数: x0name,x0init,x0ifit 确定: xfit,xfit_index,xfit_iscale
######################################################################################
def read_eamin(filename="eamin"):
    if not os.path.isfile(filename):
        print("=" * 50 + f"> ERROR: file [{filename}] NOT exist.")
        sys.exit()
    x0name = np.loadtxt(filename, dtype=str, skiprows=0, unpack=True, usecols=(0))  # 第1列: 所有势参数的名称
    x0init, x0ifit = np.loadtxt(filename, dtype=float, skiprows=0, unpack=True,
                                usecols=(1, 2))  # 第2/3列: 所有势参数初始值/势参数是否拟合
    xfit = []  # 拟合参数列表
    xfit_index = []  # 拟合参数在x0中的顺序索引
    xfit_iscale = []  # 是否归一化xscale
    for i in range(len(x0init)):
        if x0ifit[i] > 0.9:
            xfit_index.append(i)
            if abs(x0init[i]) > 0.1 and xscale == 1:  # x归一化处理
                xfit.append(1.0)
                xfit_iscale.append(1)
            else:
                xfit.append(x0init[i])
                xfit_iscale.append(0)
    xfit = np.array(xfit)  # 拟合参数列表
    return x0name, x0init, x0ifit, xfit, xfit_index, xfit_iscale


######################################################################################
# 定义Logger(), 设置标准输出到屏幕和out.txt
######################################################################################
class Logger(object):
    def __init__(self, filename="out.txt"):
        self.terminal = sys.stdout
        self.log = open(filename, "w", encoding="utf-8")

    def write(self, message):
        self.terminal.write(message)
        self.log.write(message)
        self.log.flush()

    def flush(self):
        pass


######################################################################################
# 势函数拟合程序
######################################################################################
def pot_fit(xfit):
    global ncall
    ncall = 1
    t0 = time.perf_counter()
    
    # 初始化output和ysum变量
    output = None
    ysum = 0.0
    ##########################################################################
    # 开始拟合, 打印输入文件到屏幕 (for backup)
    ##########################################################################
    print(
        f"Options: --method={method} --dstep={dstep} --dmax={dmax} --acc={acc} --maxfun={maxfun} --xscale={xscale} --yscale={yscale} --nproc={nproc} --nebimage={nebimage}")
    print("#" * 105)
    print("# Starting fitting ......")
    print("#" * 105)
    print(f"LAMMPS_0: {lammps[0]}")
    print(f"LAMMPS_1: {lammps[1]}")
    print(f"LAMMPS_2: {lammps[2]}\n")
    print(f"{' Backup: eamin ':#^105s}\n" + open("eamin", "r", encoding="utf-8").read())
    print(f"{' Backup: ref ':#^105s}\n" + open("ref", "r", encoding="utf-8").read())

    ##########################################################################
    # 设置拟合参数的上下边界
    ##########################################################################
    lower_boundary = xfit - dmax
    upper_boundary = xfit + dmax
    x_bounds = [[xfit[i] - dmax, xfit[i] + dmax] for i in range(len(xfit))]

    ##########################################################################
    # method (1) least_squares: 求解带有变量边界约束的非线性最小二乘问题
    ##########################################################################
    if method == 1:
        output = optimize.least_squares(calfun, xfit, bounds=(lower_boundary, upper_boundary), diff_step=dstep,
                                        ftol=acc, xtol=acc, gtol=acc, max_nfev=maxfun)
        ysum = calfunsum(output.x)

    ##########################################################################
    # method (2) leastsq: 最小化一组方程的平方和
    ##########################################################################
    elif method == 2:
        output = optimize.leastsq(calfun, xfit, full_output=False, col_deriv=True, ftol=acc, xtol=acc, gtol=acc, maxfev=maxfun,
                                  epsfcn=dstep, factor=100)
        ysum = calfunsum(output[0])

    ##########################################################################
    # method (3) minimize: 一个或多个变量的标量函数的最小化
    ##########################################################################
    elif method == 3:
        output = optimize.minimize(calfunsum, xfit, bounds=x_bounds, tol=acc)
        ysum = calfunsum(output.x)

    ##########################################################################
    # method (4) differential_evolution: 差分进化寻找多元函数的全局最小值
    ##########################################################################
    elif method == 4:
        output = optimize.differential_evolution(calfun, bounds=x_bounds, maxiter=maxfun, tol=acc)
        ysum = calfunsum(output.x)

    print(f"Final fitting data:\n{output}\n")
    print(f"Fitting FINISHED.  Total CPU Time: {time.perf_counter() - t0:.2f}s, ysum = {ysum:.12e}")
    return


####################################################################################################################
# calfun: 势函数计算性质
####################################################################################################################
def calfunsum(xfit):
    ylist, ysum = calfunction(xfit)
    return ysum


def calfun(xfit):
    ylist, ysum = calfunction(xfit)
    return ylist


def calfunction(xfit):
    global ncall
    t0 = time.perf_counter()

    ##########################################################################
    # 根据xfit生成x0数组, 归一化处理
    ##########################################################################
    x0 = copy.deepcopy(x0init)
    for i in range(len(xfit)):
        j = xfit_index[i]
        if xfit_iscale[i] > 0.5:  x0[j] = xfit[i] * x0init[j]  # xscale=1
        if xfit_iscale[i] < 0.5:  x0[j] = xfit[i]  # xscale=0

    ##########################################################################
    # 输出所有势参数到eam文件
    ##########################################################################
    fileeam = open("eam", "w", encoding="utf-8")
    for i in range(len(x0)):
        print(f"{x0name[i]:<15s}  {x0[i]:20.12e}  {x0ifit[i]:5.2f}", file=fileeam)
    fileeam.close()

    ##########################################################################
    # 输出优化的势参数表格
    ##########################################################################
    print("=" * 105)
    print(
        f"| {"nx":<3s}| {"Name":<15s} | {"Initial x":^21s} | {"Fitting x":^21s} | {"Difference x":^20s} | {"xscale":^6s} |")
    print("=" * 105)
    for i in range(len(xfit)):
        j = xfit_index[i]
        print(
            f"| {i + 1:<3d}| {x0name[j]:<15s} | {x0init[j]:21.12f} | {x0[j]:21.12f} | {x0[j] - x0init[j]:20.12e} | {xfit_iscale[i]:^6.1f} |")
    print("=" * 105)

    ##########################################################################
    # 生成新的势文件
    ##########################################################################
    if 'eamgen' in globals():
        eamgen(x0)
    else:
        print("警告: eamgen函数未定义，无法生成势文件")

    ##########################################################################
    # 调用LAMMPS计算, 返回给ysum, 输出计算结果表格
    ##########################################################################
    ylist, ysum = run_lammps_inFiles(ncall, time.perf_counter() - t0)

    ncall += 1
    return ylist, ysum


######################################################################################
# 调用lammps计算inFiles下的所有in文件, 处理数据, 输出data文件, ndigits为计算结果保留位数
######################################################################################
def run_lammps_inFiles(ncall=1, tgen=0.0, ndigits=8):
    t0 = time.perf_counter()
    calc_dict = copy.deepcopy(Preset_dict)
    os.chdir("inFiles")

    ##########################################################################
    # 合并in*文件为一个文件, LAMMPS计算效率会高很多
    ##########################################################################
    infile_combine = "ALL-in-combine.lmp"
    if ncall == 1:  # 只在首次运行
        file = open(infile_combine, "w", encoding="utf-8")
        for infile in os.listdir():
            if infile[0:2] == "in":
                lines = open(infile, "r", encoding="utf-8").readlines()
                file.write("#" * 105 + f"\n# {infile}\n" + "#" * 105 + "\n")
                file.writelines(lines)
                file.write("\nclear\n\n")
        file.close()

    ##########################################################################
    # 调用LAMMPS计算, 对in*文件循环计算
    ##########################################################################
    for infile in os.listdir():
        if infile[0:len(infile_combine)] == infile_combine:
            loglines, msg = run_single_lammps(infile, lamver=0)
            ###################### 分析in*结果, 截取log文件fitting结果
            for line in loglines:
                if "===fitting" in line:
                    words = line.split()
                    prop = words[0]
                    if prop not in calc_dict.keys():
                        calc_dict[prop] = round(float(words[1]), ndigits)
                    else:
                        print("=" * 50 + f"> ERROR: Repeat Calculation data: {prop} in {infile}")
                        sys.exit()
            print(f"{msg}: Number of Calculation data = {len(calc_dict) - len(Preset_dict)}")

    ##########################################################################
    # RDM计算, 对RDM*文件循环计算
    ##########################################################################
    t1 = time.perf_counter()
    for infile in os.listdir():
        if infile[0:3] == "RDM":
            loglines, msg = run_single_lammps(infile, lamver=lammps[3])
            ####### 分析RDM结果, 截取log文件fitting结果
            rdm = []
            words = []
            for line in loglines:
                if "===fitting" in line:
                    words = line.split()
                    rdm.append(float(words[1]))
            if len(rdm) > 0 and len(words) > 0:
                prop = words[0]
                if prop not in calc_dict.keys():
                    calc_dict[prop] = round(min(rdm), ndigits)
                else:
                    print("=" * 50 + f"> ERROR: Repeat Calculation data: {prop} in {infile}")
                    sys.exit()
                print(f"{msg}: Min={min(rdm):8.3f}, Max={max(rdm):8.3f}, (n={len(rdm)})")

    ##########################################################################
    # NEB计算, 对NEB*文件循环计算, 依赖于mpirun
    ##########################################################################
    t2 = time.perf_counter()
    for infile in os.listdir():
        if infile[0:3] == "NEB" and lammps[4] == 1:  # 如果不支持NEB计算, 跳过相关计算
            loglines, msg = run_single_lammps(infile, lamver=2)
            ###################### 分析NEB结果
            words = loglines[-1].split()
            PEN = [float(words[i]) for i in range(10, len(words), 2)]
            Em = max(PEN) - min(PEN)
            prop = f"{infile}"
            if prop not in calc_dict.keys():
                calc_dict[prop] = round(Em, ndigits)
            else:
                print("=" * 50 + f"> ERROR: Repeat Calculation data: {prop} in {infile}")
                sys.exit()
            print(f"{msg}: Em ={Em:8.3f}, (nebimage={len(PEN)})")

    ##########################################################################
    # Eb计算, 从ref_Eb读入计算公式: "+-"=A+B-C  "-"=A-B 
    ##########################################################################
    t3 = time.perf_counter()
    for words in ref_Eb:
        prop = words[0]
        if prop not in calc_dict.keys():
            if words[3] == "+-" and words[4] in calc_dict.keys() and words[5] in calc_dict.keys() and words[
                6] in calc_dict.keys():
                calc_dict[prop] = round(calc_dict[words[4]] + calc_dict[words[5]] - calc_dict[words[6]], ndigits)
            elif words[3] == "-" and words[4] in calc_dict.keys() and words[5] in calc_dict.keys():
                calc_dict[prop] = round(calc_dict[words[4]] - calc_dict[words[5]], ndigits)
        else:
            print("=" * 50 + f"> ERROR: Repeat Calculation data: {prop} in ref_Eb")
            sys.exit()

    ##########################################################################
    # 以表格的形式输出拟合性质计算结果
    ##########################################################################
    i = 0
    ylist = []
    fdata = []
    print("=" * 105)
    print(
        f"| {"M":<3s}| {"Properties":<15s} | {"Calculation":^12s} | {"Reference":^10s} | {"Diff":^12s} | {"Diff_pct":^8s} | {"Weight":^8s} | y (yscale={yscale:d}) |")
    print("=" * 105)
    for prop in calc_dict.keys():
        if prop in ref_dict.keys() and ref_dict[prop][1] > 0:
            i += 1
            ycal = calc_dict[prop]  # LAMMPS计算值
            [yref, weight] = ref_dict[prop]  # 参考值
            ydiff = ycal - yref  # 差值
            if abs(yref) <= 0.05:
                yratio = 500 * (ycal - yref)  # 差值比例(yscale=1): 参考值<0.05时, 给出差值*500, 否则yratio可能巨大
            else:
                yratio = 100 * (ycal - yref) / yref
            if ycal * yref < 0:   weight *= 100  # 如果ycal/yref正负号相反，权重加百倍 (有时候有用)
            if yscale == 1:
                yw = yratio * np.sqrt(weight)  # 返回差值比例, 归一化处理, 权重开平方
            else:  # yscale == 0
                yw = ydiff * np.sqrt(weight)  # 返回差值, 权重开平方
            ylist.append(yw)
            print(
                f"| {i:<3d}| {prop:<15s} | {ycal:12.4f} | {yref:10.4f} | {ydiff:12.4f} | {yratio:8.1f} | {weight:8.1f} | {yw:12.2f} |")
            fdata.append(
                f"{prop:<15s}  {ycal:20.12f}  {yref:20.12f}  {weight:20.12f}  {ydiff:20.12f}  {yratio:20.12f}\n")
        elif prop not in ref_dict.keys():
            print("=" * 50 + f"> Warning: {prop} NOT found in Reference")

    os.chdir("..")
    open("data.txt", "w", encoding="utf-8").writelines(fdata)
    ylist = np.array(ylist)
    ysum = np.sum(ylist ** 2)
    t4 = time.perf_counter()
    print("=" * 105)
    print(
        f"CPU Time: (0) eamgen: {tgen:<7.2f} (1) IN: {t1 - t0:<7.2f} (2) RDM: {t2 - t1:<7.2f} (3) NEB: {t3 - t2:<7.2f} (4) Total: {t4 - t0 + tgen:<7.2f}")
    print(f"Callfun Number: {ncall:<7d}   METHOD: {method}   ysum: {ysum:.12e}\n\n")
    return ylist, ysum


def run_single_lammps(infile, lamver=0):
    t0 = time.perf_counter()
    if not os.path.isdir("../tmp_log"): os.mkdir("../tmp_log")
    logfile = f"../tmp_log/log-{infile}"
    screen = f"../tmp_log/screen-{infile}"

    ##########################################################################
    # 调用lammps计算: 有时读取不到EAM势文件, 会出错, 需要多试几次
    ##########################################################################
    for ntry in range(10):
        ierror = os.system(f"{lammps[lamver]} -in {infile}  -echo screen  -screen {screen} -log {logfile}")
        if not ierror: break  # 没有错误就跳出循环，否则暂停一会重试
        print("=" * 50 + f"> ERROR: ({ntry + 1}) run {infile}")
        time.sleep(1.0)

    ##########################################################################
    # 确定log文件列表loglist: RDM(lamver==1)并行计算时，有多个log文件
    ##########################################################################
    loglist = [logfile]  # 默认值
    if lamver == 1: loglist = [f"{logfile}.{i}" for i in range(nproc)]

    ##########################################################################
    # 读取log文件: 尝试多次读取防止文件被占用
    ##########################################################################
    loglines = []
    for log in loglist:
        if os.path.isfile(log):
            lines = []  # 初始化lines变量
            for ntry in range(20):
                try:
                    with open(log, "r", encoding="utf-8") as file:
                        lines = file.readlines()
                    break
                except Exception as e:
                    print(f"=" * 50 + f"> ERROR: ({ntry + 1}) read {log}: {str(e)}")
                    time.sleep(0.5)
                    continue
            # 只有在lines被成功赋值且不为空时才处理
            if lines:
                for line in lines:
                    if "ERROR" in line:
                        print("=" * 50 + f"> ERROR: in {log} {line.strip()}")
                loglines.extend(lines)
    msg = f"LAMMPS_{lamver}: {time.perf_counter() - t0:7.2f}s @ {infile:<30s} "
    return loglines, msg


if __name__ == "__main__":
    sys.stdout = Logger("out.txt")
    method, dstep, dmax, acc, maxfun, xscale, yscale, nproc, nebimage = get_fit_options()
    ref_dict, ref_Eb, Preset_dict = read_ref()
    lammps = get_lammps_path()

    ###################### 调用拟合
    if method > 0:
        x0name, x0init, x0ifit, xfit, xfit_index, xfit_iscale = read_eamin("eamin")
        # 确保eamgen已定义
        if 'eamgen' in globals():
            pot_fit(xfit)
        else:
            print("错误: eamgen函数未定义，无法继续拟合")
    ###################### 根据现有的势文件直接计算
    else:
        run_lammps_inFiles()
